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Abstract 

A phase-field model for the Hele-Shaw flow of non-Newtonian fluids is developed. It extends a previous 
model for Newtonian fluids to a wide range of shear-dependent fluids. The model is applied to perform 
simulations of viscous fingering in shear- thinning fluids, and it is found to be capable of describing the 
complete crossover from the Newtonian regime at low shear rate to the strongly shear-thinning regime at 
high shear rate. The width selection of a single steady-state finger is studied in detail for a 2-plateaux shear- 
thinning law (Carreau law) in both its weakly and strongly shear-thinning limits, and the results are related 
to previous analyses. In the strongly shear-thinning regime a rescaling is found for power-law (Ostwald- 
de-Waehle) fluids that allows for a direct comparison between simulations and experiments without any 
adjustable parameters, and good agreement is obtained. 

PACS numbers: 
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I. INTRODUCTION 



The Saffman-Taylor instability occurs when a fluid is pushed by another one of lower viscos- 
ity in a confined geometry, such as porous media or a Hele-Shaw cell. It leads to the emergence 
of complex interfacial patterns whose shape is reminiscent of fingers. The study of this phe- 
nomenon has helped to establish much of our current knowledge on the self-organization 
of branched patterns 29, 3(3]. Indeed, viscous fingering can be studied under well-controlled 
conditions in the laboratory using Hele-Shaw cells, where the flow is confined to a narrow gap 
between two parallel plates. In this geometry, the full flow can be well descibed by an effec- 
tive two-dimensional problem, which greatly simplifies both theoretical analysis and numerical 
simulations. 

For two Newtonian fluids of strongly different viscosities, our understanding is fairly complete. 
In a channel geometry, the instability of a flat interface and the subsequent evolution results in the 
formation of a single finger, the so-called Saffman-Taylor finger 09I1 . Its relative width with 
respect to the channel is selected by a subtle interplay between viscous dissipation and the surface 
tension of the interface, which acts as a singular perturbation. In a radial geometry, where the 
low-viscosity fluid is injected through a central inlet, fingers are not stable and exhibit repeated 
tip-splitting to form highly ramified patterns [36]. 

Much less is known about viscous fingering in non-Newtonian fluids. Numerous experiments 
have revealed that a wide variety of patterns can be formed, including finger patterns close to the 
ones found in Newtonian fluids with either narrowing or widening of the fingers, straight fingers 
in a radial geometry that do not exhibit tip splitting, and patterns that form angular branches and 
sharp tips, reminiscent of crack networks (for a review, see y4Q). 

It is clear that the selection of these patterns is governed by the nonlinearity of the fluid itself. 
More precisely, there is a complex interplay between the geometry of the finger, which determines 
the local flow pattern. The latter, in turn, modifies the properties of the fluids. In the particular 
case of shear-thinning fluids, the dependency of the fluid viscosity on the local shear rate, which 
strongly varies in the vicinity of a finger tip, can create an effect which is akin to an interfacial 



anisotropy. The latter is known both from experiments 



fifl 



H and theory H5|,|28|] to profoundly 



affect pattern selection. Its presence suppresses tip-splitting and favors the emergence of dendritic 
patterns with sidebranches. The transition from branching fingers to dendrites observed in liquid 

o nr 

crystals [10] can thus be explained, at least qualitatively [19, 



220. Furthermore, it is not surprising 
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to see crack-like patterns in viscoelastic fluids Bill , since a high shear around the tip pushes the 
fluid into the elastic regime. 

For a more detailed and quantitative investigation of this relation between morphologies and 
the rheological properties of non-Newtonian fluids, precise numerical models would be very help- 



ful. However, in mathematical terms, viscous fingering is norma 
problem, which is quite difficult to handle numerically 11171 | l8 



ly formulated as a free boundary 



25, 



4011 . To our best knowledge, 



simulations of non-Newtonian viscous fingering using such methods have remained limited to the 
case of shear-thinning fluids in the weakly shear-thinning limit [|17l ll8ll. To overcome the difficul- 
ties due to moving interfaces, diffuse-interface and phase-field methods have become popular in 



many different fields 



ii 



23 



2611 . In phase-field models, a continuous scalar field, the phase 



field, is introduced to distinguish between the two domains occupied by the two fluids. All prop- 
erties of the fluids are interpolated through the diffuse interface, and the motion of the phase field 
is coupled to the equations of fluid dynamics. The original free boundary problem is obtained in 
the limit of vanishing interface thickness. While this approach introduces an additional scale (the 
interface thickness) into the problem, it removes the difficulties due to explicit interface tracking 
(non-uniform length change of the interfqce, topological changes). Therefore, its implementation 
is straightforward. 

In this paper, we develop a phase-field model for Hele-Shaw flow in a wide class of fluids with 
a shear-dependent viscosity, by combining a phase-field model for Newtonian viscous fingering 
previously developed by one of us [2(1 2l\ with a rigorous procedure for obtaining a generalized 
Darcy's law for non-Newtonian fluids developed by Fast et al. [1170. The model is implemented 
using a finite-difference scheme in conjunction with a standard SOR solver for the pressure equa- 
tion. We validate our model and implementation by a detailed comparison of the Newtonian case 
to the known sharp-interface solution. This allows us to estimate the errors that are due to the finite 
interface thickness and the discretization. 

Although our model is capable of describing two non-Newtonian fluids with general shear- 
dependent viscosity laws, we limit ourselves to shear-thinning fluids pushed by a Newtonian 
fluid. Indeed, this is the setting where the most precise knowledge on pattern selection in non- 
Newtonian fluids is already available, and therefore it constitutes an excellent testing ground for 
our model. Data on the shape and width of steady-state fingers for shear-thinning fluids with a 



well-characterized viscosity law have been published 1132 . 



3311 . Furthermore, these data are in 



good agreement with theoretical studies that predict a narrowing of the steady- state fingers with 
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respect to the Newtonian case [2, 



37] 



We perform simulations for two different viscosity laws, namely, a two-plateau law used in the 



simulations of Refs. 11171 . 



18], and the one-plateau law which describes well, for the experimental 



flow regime, the fluids used in experiments of Refs. B21 . 13311 . We study the effect of the shear 
thinning on the selection of the finger width, and demonstrate that our model is able to cover the 
complete crossover from Newtonian behavior at low speed to strong shear-thinning at high speeds. 
More precisely, the selection of the finger width can be understood in terms of two dimensionless 
parameters: the Weissenberg number We, which characterizes the strength of the shear-thinning 
effect, and a dimensionless surface tension T. In general, the finger width depends on both param- 
eters. However, it turns out that in the regime covered by the experiments 132j,|33j], the viscosity 
law can be well described by a simple power law. In this case, the finger width depends only on 
a single parameter, which is a function of We, T and the exponent of the viscosity law. In this 



i3 



3311 . which 



regime, our simulations are in good agreement with the experimental data of Refs. [3 
demonstrates the capability of our model to yield quantitatively accurate results. 

The remainder of the paper is organized as follows: Section HH presents the theoretical frame- 
work, the model and briefly discuss its numerical implementation. Results are then presented in 
Sec. Unl followed by conclusions and perspectives in Sec. [IV] 



H. MODEL 



A. Sharp-interface equations 

We consider two incompressible, immiscible fluids (labeled 1 and 2) in a Hele-Shaw cell of 
width W (x-direction), length L (y-direction) and gap b (^-direction, b <^ W < L). The less 
viscous fluid 2 is injected at one end of the cell with a fixed flow rate Q, causing outflow of fluid 
1 at the other end of the cell with a velocity C/qo = Q/(bW). The interface between the two 
fluids has a positive surface tension a. Both fluids, may have a non-Newtonian shear viscosity that 
depends on the local shear rate 7, ^(7*7), where r< is a characteristic relaxation time of fluid i; we 
furthermore suppose that both viscosity laws have well-defined Newtonian limits when 7 — > 0, 
which we will denote by 

As usual in a Hele-Shaw cell at low velocities (where inertia can be neglected), the scale separa- 
tion between the gap and the channel width makes it possible to simplify the full three-dimensional 
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flow problem by a long-wave approximation. The resulting two-dimensional problem is stated, for 
each fluid, in terms of the pressure field pi (which is constant across the gap) and the gap-averaged 
in-plane velocity «j. These two-dimensional velocity fields remain incompressible, 



V-Ui = 0, i = l,2. (1) 

Furthermore, for Newtonian fluids, the local averaged velocity is proportional to the local in-plane 
pressure gradient, a relationship known as Darcy's law. For non-Newtonian fluids, the relationship 
between and Vpi becomes non-linear, but can formally still be written as a generalized Darcy's 
law, 

b 2 

Ui = - Vpi, i = 1, 2 (2) 

12AifMVft|M°) 

where fxf 1 is an effective viscosity, which can be related to the original shear-dependent viscosity 
Ujjfij) for a large class of non-Newtonian fluids following the procedure developed by Fast et al. 



QJ/ZU , which is summarized and presented using the notations of the present work in Appendix lAl 
We have included the constants b, r» and //? in the argument of the effective viscosity to emphasize 
that this argument is indeed a dimensionless shear. The characteristic local shear rate can be 
estimated by the ratio of the gap-averaged velocity and the cell gap b; the order of magnitude 
of the velocity, in turn, is given by \Vp\/fJ^ (see Appendix lAl for details). Note that we have 
chosen to express the viscosity as a function of the pressure gradient (and not of the velocity as in 



L2, 



32, 



33, 



370) in order to formulate the model in terms of the interface geometry and the pressure 
field only. One should note that for a vanishing shear rate (|Vpi| — > 0), we have jjf — > fj%, and 
Eq. (O reduces to the standard Darcy's law. 

Since we are considering two fluid regions separated by an interface, we have to specify the 
boundary conditions at the interface: 

P2-P1 = CTK, (3) 

r-ui = f -u 2 = v n , (4) 

where k is the interface curvature (in the plane of flow), o is the surface tension and f is the unit 
vector normal to the interface pointing into fluid 1 . Equation © is simply the Laplace law, where 
the curvature of the meniscus between the plates has been omitted under the assumption that it is 
constant. Eq. |4] simply assures the impenetrability of the two fluids. 



In order to make this formulation more directly amenable to the construction of a phase-field 
model, we rewrite the above equations in terms of a single set of fields and material properties 

y, 



P = XiPi + X2P2, 
u = X1U1 + X2U2, 

eff , &7jjYp| . 

^Cff - Xl/il { Q . 

Hi 



v J; cff r ^2|Vpk 



(5) 
(6) 
(7) 



where and X2(x) are the characteristic functions of the domains occupied by the two fluids 

(that is, Xi(%) = 1 if the point x is occupied by fluid i, and otherwise). We thus reduce Eqs. d2l 
|3]and[T]) to just two: 



u 



12/Xeff 

V-u = 0, 



(8) 
(9) 



where 5s is a surface delta function (that is, a Dirac delta function located on the sharp interface 
S separating the two fluid domains [lj]). Now all fields, material properties and equations must 
be understood in the sense of mathematical distributions. As such, these equations, apart from 
their obvious limits at each side of the interface, are to be understood when integrated across 
the interface. In particular, integrating the normal projection of the velocity times the effective 
viscosity in Eq. © across the interface gives the Laplace pressure drop of Equation ©. Similarly, 
the condition of zero divergence of Equation © relates the normal and tangential components of 
the fluid velocity. The condition of incompressibility, when applied on the very interface, translates 
into impenetrability of the two fluids. 



B. Phase-field model 



In this section, we present the phase field approach to this problem. We first give a brief de- 
scription of the phase field (denoted by </>) and show how using it instead of the indicator functions, 
the flow equations ([8]) and © are modified. Then we present the evolution equation for the phase 
field and give a rationale for its construction. Finally, we comment briefly on how the phase-field 
model is an approximation of the sharp-interface model. 

The idea underlying the phase-field model is to introduce an additional field (0) that indicates 
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in which phase (here, in which fluid) the system is at a given space point. For the sake of simplicity 
and without any loss of generality, we consider that in fluid 1 (resp. 2) , <\> = 1, (resp. — 1). In 
addition, when crossing the interface the phase field exhibits a smooth front (kink) of finite width. 
In this general framework, the indicator functions and the 5 S function are approximated by 



Xi -> (l + 0)/2, 
X 2 -> (1 - 0)/2, 
<5 S -> |V0|/2. 



(10) 

(11) 
(12) 



Then, replacing xi and X2 by their smoothed expressions, the effective viscosity of Eq. © 
becomes 

(13) 



,i ^ 1 + ^ > „eff ^i|VpK | 1 -^ ..effAalVpl- 



/'i 



A*8 



Note that, as in Eq. (0), formally /z e ^ is a function of x because is a function of x. Darcy's law 
becomes 



u 



12/i eff (0) 



V0 



(14) 



where is the curvature of the interface computed using the standard expressions 



= V • r{(j>) and f(0) = V0/|V0|. 



(15) 



Note that k(0) is now defined in the entire space; f is the local normal to the <fi isosurface. Equation 
© for the incompressibility of the flow is not modified by the introduction of the phase field. Now, 
the flow problem is completely written in terms of the phase field. 

To complete the model, we have to introduce an evolution equation for the phase field. This 
evolution equation should have for solution a smooth interface that is advected by the flow To this 
purpose, we use the equation presented in 12011 and extended with success in [71, [80 to the case of 
vesicles: 

+ u ■ V0) = f(<t>) + w 2 V 2 (t> - w 2 K(<f))\ V0|, (16) 

with f = <p{l — (f) 2 ) the oposite of the derivative of the double well potential —<fi 2 /2 + 4 /4, 
a relaxation time, and w a small parameter that determines the width of the interface. In order to 
give a clear view of the equation, we first consider an oversimplified version of it with neither the 
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flow nor the curvature term, in a one-dimensional space. The stationary solutions of this equation 
are either the uniform solutions <p = ±1 or a front between a region where 0=1 and a region 
where = — 1: 

r 

— tanh — p. (17) 

Here, the signification of w appears clearly: it is the width of the interface. Now let us consider 
this equation (still without flow and without the curvature term) in two dimensions. Using a 
perturbation method, one can show that a weakly curved interface (radius of curvature p ^> w) 
between 0=1 and = — 1 is moving with a normal velocity proportional to 1 / p, the curvature of 
the interface. While this behaviour is expected in the case of phase transitions with non-conserved 
order parameters, here it is unphysical. In order to suppress this phenomenon, following [20] we 
add the curvature term which at dominant order is the exact opposite of the term induced by the 
Laplacian when considering a curved interface. Indeed, it can be shown that up to the third order 
inw/p, the term V 2 — n\ V0| is equal to the unidimensional Laplacian computed along the axis 
normal to the interface. Hence, using = tanh (with r the distance from the center of the 
interface), the right-hand side of Eq. (fT6l) . i.e. the driving force leading to unwanted interface 
movement, is equal to zero up to that third order. 

Finally, adding the term u ■ V0 makes the interface to be advected by the flow. Therefore, 
the dynamics of the phase field can be separated into two parts: a passive part that corresponds 
to the advection due to the flow and an active part that aims at restoring the hyperpolic tangent 
profile through the interface but does not bring any noticeable dynamics to the interface. With this 
principle in mind, it is clear that the relaxation time of the phase field must be fast enough so 
that the advection does not affect significantly the equilibrium profile. The particular choice of t<j> 
is discussed later. 

Now, that model equations have been written down, we want to stress that while the distribution 
formulation of the viscous fingering problem is just another way of writing down the same sharp- 
interface equations, the phase-field model is only an approximation to them. To be more specific, 
the phase field model introduces an additional length scale w, the interface thickness, which is a 
model parameter supposed to be small. To understand its meaning and the relationship between 
the phase-field approach and the sharp interface model, one can use the technique of matched 
asymptotics. Different asymptotic expansions of the phase field equations in powers of w valid in 
the bulk phases and through the interface, respectively, are written down. Then, matching them 
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order by order, at dominant order in w the original sharp interface-problem is retrieved, which 
indicates that the results of the phase-field model converge toward the solution of the original 
problem when w — > (the so called sharp-interface limit). In other words, the model is at least 
asymptotically correct. 

However, in numerical simulations, the value of w should be significantly larger than the space 
discetization and must remain finite. Therefore, to be able to retrieve quantitatively correct results, 
one needs to control the spurious effects introduced by the finite interface thickness and the con- 
vergence of the model toward the sharp-interface limit. This can be done by considering the next 



order in w in the matching procedure III, 



27|] . Then, new w-dependent terms are added to 



the sharp interface equations (this next order in the expansion is called the thin-interface limit). 
They actually signal the departure from the w —> limit and are the effect of the presence of the 
extra length scale. Physically, one expects their importance to depend on the ratio of w to the 
smallest genuine length scale present in the original sharp-interface model. This hypothesis can 



then be checked by simulations with decreasig values of that ratio 



21, 



22D. 



Here, we have written our mode 



so that, in the case of Newtonian fluids, it is mathematically 



equivalent to the one presented in Q20[|- The reason for this is that contrary to other phase field 
models JQ, Q] for viscous fingering, the asymptotic expansions of this model have been estab- 
lished Il20n and the numerical convergence has been checked by considering situations where the 
sharp interface solution is well known j2lb . Therefore, we are confident that unexpected finite 
interface thickness effects could only arise in our simulations in conjunction with the new feature 
here: the non-Newtonian character of the more viscous fluid. 



C. Dimensionless equations 

In order to nondimensionalize our equations, we first look at the relevant physical scales present 
in the flow, and then use the same scales to nondimensionalize the phase-field equation. Non- 
dimensionalized quantities will be denoted by a tilde. In a first step, we define dimensionless 
effective viscosity functions Ji'f by dividing the effective viscosity laws of the two fluids by their 
zero-shear limit values 

z*f (MVPI//4 1 ) = rtPf (bTi\v P \/rf). (is) 
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Next, since in a phase-field model there is a generalized effective viscosity valid throughout the 
system [Eq. (fT3l) 1 which interpolates between the effective viscosities of each fluid, we need to 
choose a single viscosity scale. This choice has to be adapted to the physical situation that is 
investigated. Here, we are mainly interested in the setting used in most experiments, where the 
more viscous fluid 1 is a shear-thinning liquid and the less viscous fluid 2 is air, that is, a Newtonian 
fluid of very low viscosity. Therefore, in the following we will nondimensionalize the effective 
viscosity by the zero-shear viscosity of fluid 1, //°. Since fluid 2 is Newtonian, we have fxf 1 = 
With the above choices, the nondimensionalized effective viscosity function becomes 



Veff{4>) 



HI 



1 + 0-eff 



1 - 



(19) 



where v is the ratio of the two zero-shear viscosities, 



(20) 



This ratio can be simply related to the quantity 



c = 



/'i 



i4 



(21) 



the so-called viscosity contrast (at zero shear), also widely used in the literature Q2C 

We furthermore measure velocity in units of the outflow velocity and lengths in units of the 
channel width W. The natural scale for the pressure gradient that arises from the Newtonian limit 
of Darcy's law is Yl^^J^jb 2 . This yields the new dimensionless quantities 



x 



y^yW V 



W 



V 



u 



uU n 



«(0) 



w 



b 2 



■Vp 



~w 
t — 



(22) 
(23) 
(24) 



Under this change of variables, the arguments of the dimensionless effective viscosity function 
given by Eqs. (1181191) become 



P>eB(<l>) 



1 + 



/if(We|Vp|) + 



1 



-V 



(25) 
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where the Weissenberg number We is defined by 



We 



12nU 



(26) 



In the remainder of this paper (except for Appendix [Aj), we will work in these new dimensionless 
variables and drop the tildes for simplicity. 

The incompressibility condition remains formally the same, and the dimensionless version of 
Darcy's law reads 



u 



1 

/x eff (0,We|Vp|) 



Vp + 1X0)- 



where 



b 2 a 



(27) 



(28) 



is a dimensionless surface tension. In summary, the flow equations contain three dimensionless 
parameters: the dimensionless surface tension T, the Weissenberg number We, and the viscosity 
ratio v. A more detailed discussion of these parameters and their role in the finger selection process 
is deferred to Sec. lIIFI below. 

To complete the set of dimensionless equations, we apply the same scaling to Eq. (fT6l) for the 
phase field. We obtain 



(29) 



and identify the dimensionless interface thickness e = w/W, the ratio of the interface thickness to 
the channel width. In order to reduce the number of purely computational parameters, we choose 
T0 = ew/U^. Indeed, w/U^ is the time it takes a flow of the magnitude of the base flow U^y to 
cover one interface thickness w, and the extra small e factor ensures that the phase field relaxation 
is one order in e faster than the forcing by the flow. We finally get 



e 2 d t = f^)-e 2 V 2 <P-k{<P)\V<P\-u- V0 



(30) 
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D. Incompressibility and boundary conditions 



In the simulations, Eq. (|30l) for the phase field and the fluid flow equations need to be solved 
simultaneously. The fluid flow part, in turn, implies solving Eq. (|27T ) taking into account the 
incompressibility condition, Eq. ©. There are several ways to implement incompressibility. 

One possibility is to take the curl of Eq. (|27l) . which eliminates the pressure field in the New- 
tonian case. Incompressibility is equivalent to the requirement that the flow is potential, that is, 
the velocity field can be written as derivatives of the stream function. The curl of Eq. (|271) yields a 
Poisson equation for the stream function. This strategy leads exactly to the model of Ref. [20] for 
Newtonian fluids, as desired. 

However, for non-Newtonian rheologies, the dependence of the effective viscosity on |Vp| 
implies that the pressure cannot be eliminated in this straightforward manner any more. Therefore, 
we use here a velocity-pressure formulation: We take the divergence of Eq. (1271) and use the 
incompressibility condition, which yields 



/2eff(0,We|Vp|) 



rX0)V0 



2/ieff(0,We|Vp|) 



(3D 



For a given configuration of the phase field 0, Eq. (1311) together with appropriate boundary condi- 
tions (discussed below) completely specifies the pressure field p. In the Newtoninan limit where 
the effective viscosity is pressure-independent, this equation is a Poisson equation for the pressure 
inside the interfacial regions where the phase field <p varies, and reduces to the Laplace equation 
in each bulk domain. In the non-Newtonian case, the source term is present also in the bulk, and 
an iterative Poisson solver must be used to obtain the pressure field for the given configuration of 
the phase field at each timestep. Then, the original Eq. (1271) immediately yields the velocity field 
u. This is then used in the next time step to advect the phase-field 0, as prescribed by Eq. (l30l) . 
More details about the numerical procedure are given in Appendix [B] 

Furthermore, boundary conditions for the phase and pressure fields are required at the edges 
of the channel. For simplicity, we will assume that if an interface crosses any of the boundaries, 
it will do so at a 90° angle, which implies that the derivatives of the phase-field normal to the 
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boundaries are zero (reflecting boundary conditions): 



6>,0 = (y = ±L/2) 



(32) 



d y( p = (x = ±W/2) 



(33) 



Since the lateral walls are sealed and hence u 



X 



0, we also have 



d x p = at x = ±W/2. 



(34) 



The only non-trivial boundary conditions are the pressure boundary conditions at the inlet and the 
outlet, where either the pressure or its gradient have to be prescribed. Since we have considered a 
flow with a fixed overall flow rate, we should prescribe the pressure gradient. 

At the outlet, only fluid 1 is present. If the interface remains far enough from the outlet, the 
pressure is simply a constant along the entire outlet, and the pressure gradient is directed along 
the y direction. Since, at the outlet, the dimensionless velocity is equal to (0, 1) (corresponding to 
a uniform flow with velocity along the y direction), the Darcy law (eq. [271) implies that the 
pressure gradient is the solution of the equation 



where the velocity enters the equation through the Weissenberg number. This equation can be 
solved numerically in a straightforward way. In our simulations, we start with an initial guess for 
the pressure gradient, which is then updated at each time step with the value found by the pressure 
solver in the vicinity of the outlet. This procedure rapidly converges to the fixed point which is the 
solution of Eq. d35l) . 

As for the inlet, we consider the case where both fluids are present. For a well-developed 
steady- state Saffman-Taylor finger, the sides of the finger are parallel to the channel walls up 
to a correction that decays exponentially with the distance from the finger tip. Therefore, if a 
sufficiently long portion of the finger is inside the simulation box, the interfaces that cross the inlet 
can be considered flat and normal to the boundary, and the fluid velocity along the x direction is 
zero in both fluids. Therefore, there is no pressure gradient along the x direction, which of course 
implies that the pressure gradient is directed along y, and constant along the inlet. 



\d y v\ = /^eff(0 = +l,We\d y p\). 



(35) 
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In contrast, the fluid velocity varies along the inlet, since the viscosity does change when cross- 
ing the interface. However, its integral along x, f^j^ u y dx, which represents the net inward flow, 
must be equal to unity, since the flow is incompressible and the fluid exits the outlet at a rate of 
unity in our dimensionless variables. Integrating the y component of Eq. ([27]) along the inlet, we 
thus obtain 

M = r+V* \ J x < (36) 

J-l/2 ^ff^Wela^l)"^ 

which constitutes a closed equation for the desired value of \d y p\ at the inlet. 



E. Simulation procedure 

In our numerical studies, our main focus is on steady-state fingers. Although we could start 
each simulation with a weakly perturbed flat interface and let it follow its natural dynamics until a 
steady finger stabilizes, this is not the most efficient procedure for parametric studies of the finger 
width as a function of T and We. Therefore, we instead first calculated an initial finger profile 
for values of the control parameters where convergence can be easily achieved, and then use the 
resulting steady-state pressure and phase fields as initial condition for a run with slightly different 
control parameters. Increasing or decreasing T and/or We in small steps, we are thus able to follow 
the steady-state solution branches over a substantial parameter range. 

When performing the first computation for a given viscosity law, we set the initial interface pro- 
file to a semi-elliptic bubble (of width W/2 and length W) growing from the inlet of the channel. 
The initial configuration of the phase field is a hyperbolic tangent profile in the elliptic coordinates, 
and its zero contour is located at the elliptic bubble interface. The simulations are performed in a 
channel with a length of L = 5W. The bubble increases in size and depelops into an elongated 
finger. When it reaches a reference position (typically, located at twice the channel width from 
the inlet), the whole domain is translated backward by one grid spacing (in other words, the finger 
is pulled back by one grid point). The velocity of the finger is computed by measuring the time 
between two successive pullbacks. The finger width is measured at the entrance of the channel 
when a pullback occurs. We consider the steady reached when both tip velocity and finger width 
vary less than a fixed value (here chosen to be 1CT 8 , to be compared with a typical tip velocity of 
2 and a typical finger width of 0.5) between two pullbacks. 

Values of T of the order of 10~ 2 yield a rapid convergence to a steady-state finger, both for 
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Newtonian and non-Newtonian fluids. For the latter, the convergence is more difficult to obtain 
because of the nonlinearities in the viscosity laws. Typically, we calculate the first finger with 
a low value of the Weissenberg number We for which these nonlinearities are small; We is then 
increased progressively up to the desired value. In this way, values up to 10 2 can be treated, for 
which the pressure solver would have otherwise not converged. As for T, the values we are able to 
attain are limited both from below and from above. For small values of T, results become sensitive 
to the discretization and the interface thickness, as will be detailed below. For large values of V, 
the finger width becomes close to unity, and the tail of the phase-field profile starts to interact with 
the sidewalls. 

The solutions found in our simulations are single fingers propagating at constant velocity along 
the channel. We consider fingers symmetric with respect to the channel mid-line x = 0. This 
allows us to reduce the computation time by limiting the numerical domain to half the channel: 
< x < 1/2). The validity of this procedure was checked by occasinally performing computa- 
tions in the full domain: fingers started with an axis of symmetry shifted away from the mid-line 
always relax towards the center of the channel in finite time. We have also checked that increasing 
the length L of our simulation domain (changing the aspect ratio L/W) does not change the re- 
sults. Indeed, in our typical steady-state configuration, the back of the finger is cut off at twice the 
channel width behind the tip, where its flanks are almost flat and fluid 1 is almost at rest. Further- 
more, the pressure field becomes almost linear far ahead of the tip, and a distance of three times 
the channel width is enough to resolve all non-trivial features of the velocity and pressure fields. 

In our simulations we let the finger extend inside the channel until the tip crosses a reference 
position along the y axis. When this happens, the time step is truncated so that the finger tip 
advances exactly to the pullback coordinate; the whole field is then pulled one grid step backward. 
The velocity of the finger is obtained by computing the average velocity between two successive 
pullbacks. The finger width is measured at the entrance of the channel when a pullback occurs. 
The stationary state is declared to be achieved when both tip velocity and finger width vary less 
than a fixed value, here chosen to be 1CT 8 . 

The computation time necessary to achieve the stationary state for a given B value can be 
significantly reduced when the run is initialized with a finger profile close enough to the converged 
state. Hence we applied the following procedure to obtain selection curves in the Newtonian and 
the shear thinning cases: 

For the first computation of the set, the phase field is initialized with a semi-elliptic bubble 
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growing from the inlet side of the channel. The small radius spans over half the width of the 
channel and the big radius is arbitrarily chosen to be the channel width in the longitudinal direction. 
The phase field obeys a hypertangent profile in elliptic coordinates. A moderate B parameter of 
10~ 2 is chosen to compute the first finger profile. The parameter B is then varied towards zero 
and towards infinity to move along the selection curve. Each computation is initialized with the 
finger profile at numerical convergence. Attainable B values are both limited in the small and large 
limits. In the former case, the interface thickness needs to be reduced, and thus the grid refined, in 
order to retain the relevant selection mechanism. In the latter case, the stationary solution can be 
destroyed when the phase field is too close to the boundaries. In the non-Newtonian case the first 
step is more difficult because of the nonlinearities in the viscosity. We do not impose the correct 
pressure condition at the outlet, but rather let it relax as time is stepped forward. 



F. Control parameters and finger selection 

The independent parameters that appear in our equations are the zero-shear viscosity ratio v 
(constant for a given pair of fluids), the Weissenberg number We, which controls the intensity 
of the shear-thinning effect, and the dimensionless surface tension T. It is noteworthy that in 
experiments performed with a single Hele-Shaw cell of fixed width and gap spacing, both We 
and 1/r increase linearly with [see Eqs. (|26l),(|28l)l, which is the only parameter that can be 
externally controlled. The full two-dimensional parameter space can hence only be explored in 
experiments by varying the channel geometry as well as U^. In contrast, in the simulations it 
is easy to vary these two parameters independently, and to determine the selected finger width. 
However, it is useful to take some additional considerations into account. 

It is known that two main ingredients determine the finger width: the dimensionless surface 
tension and the anisotropy of the interface or the medium. In shear-thinning fluids an effective 
anisotropy arises from the fact that the in-plane velocity and thus the shear are maximal at the tip, 
and decay when going to the sides of the finger. As a consequence, the viscosity and hence the 
mobility in Darcy's law vary along the interface. Thus, the strength and nature of this effective 
anisotropy are essentially controlled by the Weissenberg number and the functional form of the 
viscosity law. 

.et us now turn to the dimensionless surface tension. For Newtonain fluids, it was shown 



[|35|l that the selection of the finger width is determined by a single dimensionless parameter B, 
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which represents the ratio of stabilizing (capillary) to destabilizing (viscous) forces. The 
are proportional to the finger speed and the difference of the two viscosities (see e.g. Ref. |4O0). 
In shear-thinning fluids, the relevant viscosity is the one in the vicinity of the tip, and the correct 
definition of the parameter B is 

B = h ^—, . (37) 

12W 2 U tip \jj°fif (We|Vpti P |) - f£] 

Using the definition of T and the fact that mass conservation for an incompressibe fluid enforces 
Uoo = A {/tip for a steady-state finger of relative width A, we find the following relation between B 
andT: 

Z = ™ . (38) 

{/up <(We|Vp tlp |) - v /if (We|Vp tlp |) - v 

Ideally, we would like to explore the parameter space along lines of constant B in order to track 
only the influence of the effective anisotropy (the selection parameter of the isotropic Saffman- 
Taylor problem is then constant). However, B is difficult to control directly in our simulations: the 
effective viscosity at the tip, which is needed to calculate B, depends on the finger speed, which 
is itself the result of the selection to be investigated. Therefore, we explore the width selection by 
varying either We at fixed T, or T at fixed We, and calulate B a posteriori using the tip speed U t i p 
and pressure gradient | Vptip | extracted from the simulations. Note that this procedure is perfectly 
analogous to the one followed in experiments: the viscosity at the tip is estimated a posteriori using 



the measured finger speed 11321 . 13311 . Keeping B constant is more involved, and would require some 
iterative trial and error procedure, which is perfectly feasible but cumbersome. 

A last point that deserves brief mention is the viscosity ratio v. In the case of air pushing a 
viscous fluid, v is extremely small, so that the viscosity of the air can be neglected altogether. 
In our numerical formulation, however, it is difficult to simulate very small values of v, because 
Eq. (I3TI) then has extremely different numerical stiffness in the two bulk domains, which makes 
the convergence of the pressure solver delicate. In our simulations, we have typically used values 
of v ranging from 5 x 10~ 2 to 5 x 1CT 4 , which are large enough to guarantee a robust and efficient 
solution of Eq. (I3TI) . One could think that these are small enough to neglect v in the denominator of 
Eq. (l38l) . However, as will be seen below, for a viscosity law without lower bound (such as a power- 
law), the viscosity of the shear-thinning fluid will become comparable to or even smaller than that 
of the pushing Newtonian fluid even for v <C 1, for sufficiently high Weissenberg numbers. In the 
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latter case, the fingering instability disappears altogether. We insist that this is an entirely physical 
effect that should be experimentally observable in fluid couples of not too different viscosities. 



III. RESULTS 



A. Newtonian fluid 



In order to test our model formulation and its numerical implementation, we start by performing 
simulations in Newtonian fluids. The simulations converge without difficulty to a steady-state 
finger solution. In Fig. Q] we display a comparison between a typical finger shape extracted from 
our simulations and the analytical solution of Saffman and Taylor Q39Q , 
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(39) 



After fitting A, the agreement between the computed and analytical curve is good. There are some 
small discrepancies close to the finger tip that are to be expected, since the solution given by 
Eq. (|39l) does not contain the effect of surface tension. That the latter is correctly incorporated 
into our model is proven by the results shown in Fig. [2l where we display the selection curve 
for the finger width at fixed values of v and e as a function of the dimensionless combination of 
parameters 4Ai?7r 2 /(l — A) 2 used in the classical work of Mc Lean and Saffman [13511 . which we 
compare to. The agreement is excellent, except for very small values of T. This constitutes an 
extremely sensitive test for our model since the finger width is selected by the surface tension (via 
the selection parameter B) through a singular perturbation mechanism. 
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Figure 1: Comparison of computed interface ( T = 0.01, 1^=0.05, e=0.02, Ax=0.01) and analytical solution 
of Saffman and Taylor, A = 0.58. 
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Figure 2: Comparison of computed finger width A (y = 5 x 10~ 3 , e=0.02, Ax=0.01) and semi-analytical 
solution of McLean and Saffman. 

In Fig. |3] we replot the selection curve directly as a function of T, to make the deviations from 
the analytical prediction for small values of T most apparent. These deviations take place for 
T < 0.01; Below that value, the decrease of the finger width A with the dimensionless surface 
tension T to the predicted limit value of A = 0.5 for r — > (note that for a Newtonian fluid, B 
is just proportional to Y) is interrupted by a small "bump". Two effects limit the precision of our 
results. First, it is expected that at low values of B smaller values of e are needed to obtain properly 
resolved results. The reason is that the wavelength of the marginally stable mode of the linear 
Saffman-Taylor instability scales as ~ \[B. As in any phase-field model, the correct interface 
dynamics can only be guaranteed a priori when e remains smaller than this value (i.e., well into 
the thin-interface limit). Deviations from the sharp-interface solution are thus simply a sign of 
insufficient resolution of the relevant length scale by the phase field. The second effect is purely 
numerical: when T is decreased the surface tension effect becomes numerically small. More 
precisely, the pressure gradient accross the interface created by the Laplace pressure becomes 
smaller and smaller with respect to the global driving pressure gradient. Therefore, discretization 
errors can become significant. In particular, the anisotropy induced by the discretization on a 
regular lattice can have a strong effect on the solution. This is especially critical, since it is known 
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Figure 3: Steady-state finger width versus T, (a) for various e=0.02, 0.016, 0.01 and Ax=0.01, 0.008, 0.005, 
and (b) for fixed e=0.02 and increasing resolution of discretization, Az=0.01, 0.008333, 0.005. i/=0.05. 



that even a small amount of interfacial anisotropy dramatically modifies the selection mechanism 



L15, 



28] 



In Fig. [3l we test the importance of these two effects. Whereas a reduction in e at fixed res- 
olution (that is, constant e/Ax) reduces the height of the "bump", the change of sign in slope 
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occurs always at similar values of T. In contrast, if the mesh is refined at fixed e, the change in 
slope is shifted towards smaller values of T. This indicates that the numerical discretization error 
is the dominant effect. Since a further decrease in the grid spacing would require a much larger 
computation time, we have limited our study to the regime of intermediate values of T. 

Incidentally, an observation we find worth reporting is that of symmetrical pulsating fingers 
(i.e. time-periodic solutions with oscillating width and tip velocity), albeit in a parameter region 
where the numerical convergence is not guaranteed (T slightly below (1Q~\ These oscillations 
disappeared after further grid refinement. This is consistent with the picture [6|] according to which 
the threshold in the logarithm of the amplitude of the noise (here numerical and related to the grid) 
needed to nonlinearly destabilise a Saffman-Taylor finger decays linearly with — T~P, (3 > 0, 
(3 ~ 0.5. 

B. Shear-thinning fluids 

To study the effect of shear thinning, we first need to specify the viscosity law. As an example, 
we take a two-plateau Carreau fluid, whose viscosity obeys the equation 

M 7 } ~r = (1 + (r7) 2 ) ( - 1)/2 . (40) 

Besides the already introduced relaxation time r and zero-shear viscosity fi®, this law has an 
inifinite- shear asymptote at the value and an exponent n. It describes three regimes: two 
Newtonian plateaux at zero and infinite shear, where the viscosity is independent of the shear rate, 
and a shear-thinning region in between. The ratio of the heights of the two plateaux can be defined, 
a = fif I fx®, whereas the slope in the shear-thinning regime is determined by both n and a. 

In the following, we address two limiting cases of this general law: the weakly (a not too 
small, see below for a more precise statement) and the strongly (a — > 0) shear-thinning regimes. 
No analytic expression for the corresponding effective viscosity (to be used in Darcy's law) is 
known in either limit. 

1. Weakly shear-thinning fluids 



We first consider the weakly shear-thinning case and set n = —1 in Eq. (|4Q|) to make contact 
with Ref. II 1711 . It was shown there that the resulting law translates into an effective viscosity 
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in Darcy's law as long as a > 1/9, so for practical purposes that sets the minimal value of a 
which we mean when we refer to the "weakly" shear-thinning case. However, no closed analytical 
expression for this effective viscosity seemed possible, but the same functional dependence as the 



viscosity law Eq. (1401 with n = — 1 turned out Ill7ll to provide an excellent approximation for it 



f= l + a |WeV rf 
1 + |WeVp| 2 

where we recall that the Weissenberg number is given by We = 12t 1 [/ 00 /&. We therefore use this 
law in the remainder of this section. 

The value a = 1 corresponds to a Newtonian fluid; when a decreases, the viscosity variations 
become steeper. We recall that there are now two independent parameters that control the finger 
selection (on top of a): T, as in Newtonian fluids, and We, which measures the strength of the 
shear-thinning effect. We begin by investigating the role of We. 

Let us first illustrate the origin of the effective anisotropy effect for shear-thinning fluids by 
display maps of the local effective viscosity function £b e g, Fig.H in various flow regimes, i.e., for 
various ranges of We values. We find it clearer to begin with a description of the velocity field, 
since it relates directly to the local viscosity through the shear rate, which is proportional to the 
gap-averaged velocity. Far ahead of the finger, the local velocity is = 1 as in the outlet. The 
speed increases when the finger is approached, since the finger tip speed, U tip = U^/X ~ 2 is 
larger. Indeed, this is the maximal speed in the system. Further upstream (along the finger flanks) 
the speed of fluid 1 decreases to its limiting value, which can be computed using Eq.[36]and is of 
the order of u/X. For v — (inviscid pushing fluid) the limiting value is 0. 

With this picture in mind, the shear-thinning phenomenon upon increasing We should be 
clearer. For We < 1, we remain in the low-shear Newtonian plateau of the viscosity, which is 
hence homogeneous. As (We > 0.1), the speed at the finger tip enters the shear-thinning regime, 
so the effective viscosity exhibits a well-marked minimum there. Furthermore, it increases towards 
its Newtonian limit along the finger sides, and it also increases ahead of the finger and towards 
the outlet. This picture remains valid when We increases further, with the only difference that the 
region where the Newtonian regime is reached is sent further upstream along the finger flanks. 
Eventually, for (We > 5), a third regime appears: The fluid at the finger tip enters the high-shear- 
rate plateau of the viscosity law, so the viscosity becomes homogeneous in a growing region close 
to the tip, although it remains its absoute minimum in space. Soon the outlet is taken by this 
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homogeneous-viscosity region, since the speed there is typically just a factor 2 smaller than at the 
tip. However this is not the case of the finger flanks, where the fluid speed decreases much more 
upstream, so they remain a shear-thinning zone, provided v is small enough. This shear-thinning 
zone expands and moves upstream as We is furhter increased; if v > is kept constant, it will 
eventually reach the inlet, and if We is even increased further, the whole shear-thinning zone will 
"pass" through the inlet until the spot reaches the high-shear plateau and the viscosity becomes 
homogeneous again everywhere (but now lower). This happens at We ~ 1000 for v = 5.10 3 , 
regardless of the finger length simulated. 

In Fig. [51 we show the selected finger width as a function of We at fixed T for various values 
of a and T. For small values of We, the shear-thinning fluid is in the high viscosity plateau and 
the finger width is almost constant. When We is increased and approaches unity, the finger width 
decreases. The curve goes through a minimum, after which the finger width increases with We 
until it becomes constant again when the shear-thinning fluid enters the second plateau. 

The finger width for given a and T is larger at We >> 1 than at We << 1. This is a conse- 
quence of the relation between the control parameter T and the tip selection parameter B already 
discussed in Sec. Ill Ft the control parameter V is defined with the viscosity of the first Newtonian 
plateau. However, at high shear rates, the fluid around the tip is in the second plateau, and there- 
fore the width selection is governed by the corresponding value of the viscosity. Neglecting the 
viscosity of the Newtonian fluid (that is, setting v = 0), we obtain at high Weissenberg numbers 
the simple relation B = T/a, whereas for low We, B = T. Since a < 1, larger fingers are 
selected for high We. This argument is corroborated by the two curves for T = 0.02, a = 0.3 and 
T = 0.01, a = 0.15, which tend to the same finger width at high We (Fig. [5]). Indeed, they have 
the same value of B = 0.033 in that regime. 

A noteworthy feature of Fig. [5] is that all the curves for T = 0.01 exhibit finger widths that 
are lower than 0.5, which is the smallest value that can be achieved in Newtonian fluids. This 
narrowing is due to the effective anisotropy induced by the shear-thinning effect in the medium, 
as can be appreciated from the viscosity maps in Fig. 0] the region of lower viscosity right in front 
of the finger tip facilitates the advance of the interface in the center of the channel. It is thus not 
surprising that the lowest values of the finger width are reached for We ~ 1, where the variations 
of the viscosity close to the tip are the strongest. Furthermore, this effect increases with decreasing 
a, as can be seen by comparing the three curves obtained at T = 0.01 in Fig. [51 They coincide 
at small We values since the first Newtonian plateau is the same for all the curves. When We 
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approaches unity, the finger width decreases, with smaller a giving rise to narrower fingers. This 
is to be expected since a smaller a implies stronger variations of the viscosity with the shear rate, 
and thus a stronger effective anisotropy. At We rs 10, the curves cross. Now lower values of a 
give rise to wider fingers. This is due to the global decrease in viscosity in the shear-thinning fluid 
already discussed above, together with the weakening of the shear-thinning effect around the tip 
when the fluid enters the second Newtonian plateau. 

From the preceding discussion, it is clear that for the effective viscosity law given by Eq. @T] 
the strongest shear-thinning effect occurs for We ~ 1. Therefore, next, we fix We = 1 and study 
the selected finger width as a function of T for various values of a. As discussed previously, in 
order to display the results in a meaningful way, finger widths need to be plotted as a function of 
B, which can be calculated a posteriori using Eq. (T37T ). Figure [6] displays three selection curves 
for a = 0.9, 0.3, 0.15 and v = 5 x 10~ 3 , compared with the corresponding Newtonian curve at 
i/ = 5x 10~ 3 . The curve for a = 0.9 is very close to the Newtonian one; with decreasing values 
of a, the selected finger width decreases at fixed B, which is consistent with the picture of an 
effective anisotropy increasing with a. It should also be noted that, as for the Newtonian fluid, 
a "bump" occurs in the selection curve due to discretization effects; however, for strongly shear- 
thinning fluids there clearly exists a range of B for which the solution is not affected by numerical 
artifacts, and for which stationary fingers display a width A < 1/2, which would be impossible for 
a Newtonian fluid. 



2. Strongly shear-thinning fluid 

We now turn to the case in which the shear-thinning effect is strong, that is, the infinite-shear 
viscosity can be neglected in front of the zero-shear viscosity, fif — > or a — > 0. Then, the 
high-shear plateau disappears, and the two-plateau law of Eq. (|40l) becomes a one-plateau Carreau 
law (see e.g. |4J]) 

/i 1 (r 1 7)=/i?(l + (r 1 7) 2 ) (n - 1)/2 . (42) 
This expression has been shown to provide agood fit to the aqueous solution of the polymer 



xanthane used in the experiments of Ref. Q32 



3311 . In our simulations, we now use n = 0.5. 



Even for this simplified law, again no analytical expression exists for the corresponding ef- 
fective viscosity /i efr (We| Vp\) to use in Darcy's law. We have therefore tabulated the effective 
viscosity for our numerical calculations, following the procedure of Appendix [A] However, in 
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the limit of large shear rates, T17 ^> 1, Eq. (1421) reduces to the Ostwald-de-Waehle power-law 
viscosity, /ii ~ (ti^)"™ 1 , and it is easy to show that the effective viscosity asymptotically behaves 
as 

/ \ (n— l)/n 

/i efr (We|Vp|) ~ (We|Vp|J . (43) 

Let us start by discussing the effect of the Weissenberg number. The curve of the finger width 
versus We obtained at constant T = 0.01 is shown in Fig. |7J The onset of the shear-thinning 
regime occurs for We close to unity as in the weakly shear-thinning case. The finger width then 
reaches a minimum once the whole tip region is in the shear-thinning regime. When We is further 
increased, the width increases monotonously, without exhibiting a plateau as in Fig. |5] This is 
of course due to the fact that there is now no second plateau in the viscosity law itself either. 
The viscosity continues to exhibit a marked minimum at the tip, which implies that the effective 
anisotropy is present for any We > 1. At the same time, the viscosity everywhere in the shear- 
thinning fluid decreases with increasing We, which leads to an ever increasing value of the tip 
selection parameter B, and therefore to an increase in width. 

Ideally, in order to separate the global variation of the viscosity from the appearance of the 
effective anisotropy in the shear-thinning fluid, the finger width should be studied at fixed B in- 
stead of fixed T. This is difficult since B can only be evaluated a posteriori, as already discussed 
in Sec. Ill F[ However, a procedure can be devised that yields a clearer view: instead of varying 
We at constant T, one may also vary simultaneously We and T to keep constant the dimensionless 
surface tension defined with the viscosity at the outlet: 

r out = 5V s = (44) 

UW^^fif (We|Vpout|) /if(We|Vp out |) 

When We is varied, the new effective viscosity at the outlet is computed through the pressure 
gradient value there, which is the numerical solution of Eq. (|35l) . T is then chosen to keep r out 
constant. 

The rationale for this procedure is the following: In the fully shear-thinning regime, the whole 
tip is surrounded by a fluid region in which the effective viscosity scales as a simple power law. 
Intuitively, changing We in this regime should not alter the effective anisotropy at the tip, and the 
finger width selection should be governed by the only selection parameter left, the dimensionless 
surface tension B. This scenario would be perfectly consistent with the theoretical studies of 
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Refs. [|2J,[37J]. In the power law regime, the viscosities at the tip and at the outlet scale in the same 
way. Therefore, for an inviscid pushing fluid, carrying out simulations with r out constant should 
leave the value of B and hence the finger width constant. Indeed, it can be seen in Fig. [7] that 
the finger width for large Weissenberg numbers varies much less when keeping r out than when 
keeping T constant. The residual increase of A with We is due to the finite viscosity of the pushing 
fluid 2(y> 0). 

Let us see the relation between T, B, and r out in detail: In the power-law regime of Eq. (1431) . 
jy^ji/n Then, the viscosity of the Newtonian fluid is negelcted (that is, we set v = 0), 



Eq. (1381) becomes 

B = ^ L_ ^r. (45) 

U tipf ,f(We\Vp\ tip ) U% p 

Taking into account that, for a steady-state finger of width A, U tip = U^/X, and that A is a unique 
function of B, we obtain 

> U^ n T ~ We 1 ""!, (46) 



X(B) n 

and it is clear that B (and hence A) is fixed by the product rWe 1_n . 

Similarly, it can be seen that r out defined by Eq. (1441) scales as ~ TWe 1- ™. Therefore, for an 
inviscid fluid 2, keeping r out constant amounts to keeping B, and hence the finger width, fixed. 
The reason for the residual increase in A with increasing We but fixed r out in Fig. [7] is the finite 
viscosity ratio v. Indeed, this ratio, which is a constant independent of We, appears in the relation 
between B and T, Eq. (|38l) . Therefore, as the viscosity of the shear-thinning fluid 1 decreases with 
increasing We, the denominator gets smaller. As a result, even at fixed r out , B increases with 
increasing Wejeading to an increase in the finger width A. 

This analysis shows that the dimensionless surface tension B in the power-law regime is deter- 
mined by the parameter TWe 1_n . Our intuition that the effective anisotropy remains constant in 
that regime suggests that the whole dynamics is controlled by this single parameter. Substituting 
the expression for the power-law effective viscosity, Eq. (1431) into Darcy's law, Eq. (|27T) . we get 



u ~ - (We|Vp|) 



v p + r«(0)^ 



(47) 



Assuming that we have a solution of the problem (i.e.: velocity field, pressure field and finger 
shape, implicitly given by (f>) for a given set of We and T, we consider a situation where the 
product rWe 1_n is kept constant while We is multiplied by a positive value £ (this amounts to 
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multiply T by £ n_1 ). In this case, considering eq@7l it is clear that if Vp is also multiplied by 
£ n_1 , the velocity field will be kept unchanged (and thus obey the boundary conditions and the 
incompressibility condition for fluid 1). In the case where fluid 2 is inviscid, its pressure gradient 
is zero, so the above rescaling for the pressure gradients in fluid 1 yields indeed a strictly valid 
solution in the whole domain and the dynamics depends only on the reduced parameter rWe 1_n . If 
the viscosity of fuid 2 is finite (but negligible), < v « 1, this rescaling is a good approximation. 

These predictions are indeed borne out by the simulation results. In Fig. [81 we show the relative 
finger width as a function of the only relevant parameter rWe 1_n (here, n = 0.5) for two series 
of simulations carried out at different values of We. The two simulation curves, which differ by a 
factor of 5 in the Weissenberg number, superimpose almost perfectly. 

It turns out that the high-shear limit We| Vp| 3> 1 is also the relevant regime for the description 
of the experiments of Refs. D2l 13311 . Therefore, we show in the same plot experimental data for 
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various channel geometries. We have included complete data sets from Refs. H32I . 
exhibit first a decrease in the finger width with decreasing control parameter, but below a certain 
value they start to increase again, contrary to the theoretical predictions. This was attributed later 
tl2D to the onset of inertial effects, which are obviously not contained in our model. Nevertheless, 
we have included all data points in our plot in order to avoid an arbitrary cutoff. The part of the 
data not affected by inertia (the part with a positive slope) is quite close to our numerical curve. It 
is interesting to note that to rescale the experimental data, only the channel geometry and the flow 
rate Q (or, equivalently, the finger speed and width) have to be known; no data on the viscosity 
of the tip are needed. Furthermore, it is useful to stress that the scaling analysis makes it possible 
to meaningfully compare simulations and experiments, even though they are not carried out at the 
same parameters. The experimental data correspond to high Weissenberg numbers (We ~ 10 3 ) 
and extremely small values of the viscosity ratio v (since the pushing fluid is air); carrying out our 
simulations at these parameters would have been quite a numerical challenge. 



IV. CONCLUSIONS 



We have developed and validated a phase-field model for viscous fingering in shear-thinning 
fluids in a Hele-Shaw cell. It can be used for fluids with arbitrary shear-dependent viscosity, 
provided that the viscosity function is not too steep to allow for the calculation of the effective 
viscosity function by the method described in Appendix El We have also shown that the model 
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is capable to describe the full crossover from Newtonian to strongly shear-thinning behavior, and 
to make quantitative contact with experimental results. It is therefore a useful and robust tool for 
further investigations of the precise relationship between the rheology of the shear-thinning fluid 
and the pattern formation process. 

We have investigated the selection of the finger width in the channel geometry for two different 
shear-thinning laws. One is the model fluid already used in Ref. [17] that exhibits two plateaux in 
the viscosity function at low and high shear rates, and describes weakly shear-thinning fluids. We 
have found that a narrowing of the fingers below the limit A = 1/2 for Newtonian fluids is observed 
only in the regime where most of the variations of the viscosity occur in the vicinity of the tip. This 
confirms the idea that the self-organization of the medium provides an effective anisotropy leading 
to sharper finger tips. The second rheological law investigated describes well the behaviour of the 
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strongly shear-thinning fluids used in the experiments of Refs. 11321, 
power-law viscosity at large shear rates. In this case, the system reaches a scaling regime where 
the finger width depends on a single parameter, simply expressed in terms of the channel geometry 
and the exponent of the viscosity law. This scaling makes it possible to compare simulations and 
experiments, even though they are not carried out at the same parameters. Reasonable agreement 
is obtained. 

In the future, it would be interesting to use this model for a systematic investigation of pattern 
selection as a function of the viscosity law, especially in the regime of narrow fingers. However, to 
attain this "needle regime", improvements in the numerical algorithm will be needed, in particular 
a refinement of the grid spacing at the interface. This could be achieved using adaptive meshing 
algorithms. Finally, the model can also be used without any difficulties to simulate fingering in 
radial Hele-Shaw cells and to study the transition from tip-splitting to stable dendritic growth. 
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Appendix A: DARCY'S LAW FOR NON-NEWTONIAN FLUIDS 



A (Newtonian) viscous fluid in a Hele-Shaw cell or porous medium obeys Darcy's law: its 
velocity is proportional to the local pressure gradient for not too high gradients, since then inertia 
can be neglected. The proportionality constant can be understood as a mobility, and it depends 
on the fluid viscosity and the characteristics of the medium. In particular, in a Hele-Shaw cell 
these "medium" characteristics are purely geometrical, since the mobility appearing in Darcy's 
law is actually an average across the cell gap. The underlying idea is to project the actual three- 
dimensional problem into and effective bidimensional problem in the plane of the glass plates, 
taking advantage of the fact that the cell gap b is much smaller than any other length scale in the 
problem. In this projection procedure, one starts from the Stokes equation for any fluid labelled 
by i, 

V • {ViVu) = Vp. (Al) 

All quantities have their corresponding dimensions; some of their dimensionless counterparts, as 
defined in particular by Eq. (T23T) . will only appear at the end of this Appendix and will then be 
denoted by a tilde on top of their respective symbols. 

In the left hand side of this Stokes Equation (IA1I) . d x and d y are neglected with respect to d z , 
much stronger due to the small gap thickness, and one considers only the in-plane flow (x and y 
directions). We continue to denote the bidimensional versions by u and V to keep the notations 
simple. Note that, here, u Is a function of z. Integrating once, one gets 

Hid z u = zVp. (A2) 

Darcy's law is then obtained by integrating once more to get the in-plane velocity u and av- 
eraging the latter over the cell gap. While this is straightforward for Newtonian fluids where the 
viscosity is just a constant, in the non-Newtonian case where the viscosity depends on the shear 
\d z u\, this is only possible if this function is invertible. In the following, we detail the steps to 
obtain Darcy's law in this case, in the spirit of Fast et al. Jl7 1: 

We rewrite the viscosity as fa = $fa(rf\d z u\ 2 ), where /x? is the zero-shear viscosity and fa 
is a general, dimensionless viscosity function of a dimensionless argument, with some internal 
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relaxation time of the fluid. We take the modulus of Eq. (|A2I) and multiply it by r, to get 



fii{s 2 )s = (z, (A3) 

where s = Ti\d z u\ and £ = Xj|Vp|//i°. As long as s/Ij(s 2 ) is an invertible functional], 
this equation constitutes an implicit function s 2 (( 2 z 2 ), which we reinject into fii(s 2 ) to get 
fii(s 2 (( 2 z 2 )) = nl(( 2 z 2 ). We can now formally solve Equation (IA2I) for d z u: 

and integrate it to get the in-plane velocity 



« = — / ^TTT^T (A5) 
-6/2 



^ 7-6/2 K(C 2 ^ 



Finally, we compute the gap-averaged velocity 

1 Z^ 2 

(m) = - / wcte. (A6) 

Oj-b/2 

After performing this latter integral by parts and taking into account that the integrand is even, we 
obtain 

Vp f b ' 2 z 2 dz 



/4(C 2 * 2 )' 

At this point, we have obtained a relationship between the gap-averaged velocity and the pres- 
sure gradient which is non-linear since the pressure gradient appears not only in the prefactor, but 
also in the integral (in the form of the factor Q. This relation can then be used to define an effec- 
tive viscosity that depends on the pressure gradient. For computational purposes it is preferrable 
to change the variable of integration from z to s according to Eq. (IA3I) . In doing so, we go back 
from the inverse function /4(£ 2 z 2 ) to the original shear viscosity function /Ij(s 2 ). We get that 

1/2 «"* ■/" (( ^ (A8) 



o C 3 io " da 

K ... ,, fWi|Vp| 



where x = with b( = — l — — . (A9) 
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Integrating by parts once more we finally obtain 



(u) = - b -^(bcr 3 <U 2 (x 2 )x 3 - f x tfpyds} . (mo) 

A*i L Jo ) 

This can be formally rewritten in the form of a Darcy's law, 

w-idSPb <A11) 

where 10 _* = &C -3 (/*i 2 (xV " f* fii 2 (s 2 )s 2 ds) , (A12) 



12/if(&C) 
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with a mobility where the purely geometrical factor of 12 for Newtonian fluids has been replaced 
by a complicated function of the variable b(. This variable actually represents the dimensionless 
shear. Rewriting ( in terms of the original quantities, and then scaling the pressure as in the main 
text [Eq. (1231)1, it becomes 

bn\Vp\ UnUoo rf^ j We\fp\ if i = l 

K = — = — r o\^P\ = \ ~-> , ( A13 ) 

fH " ri I (r/v)We\Vp\ if % = 2 

where the Weissenberg number We and the zero-shear viscosity ratio v are defined by Eqs. (|26l) 
and d20l) in the main text, and 

r = —, (A14) 

n 

is the ratio of the characteristic time scales of the two fluids. Note that we have here allowed 
for two different shear-dependent viscosity laws. The global interpolated effective viscosity law 
becomes then 



VM) = -^—lA ( 1 |Vp|) + — ^— VH 2 ( ^|I CV PI) 



= l ~^Pf (We|^p|) + l_i^f (We|^|r/z/). (A15) 

The formulas of the main text, valid if fluid 2 is Newtonian, can then be obtained by setting r = 
and jlf 1 = 1 . 
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Appendix B: NUMERICAL METHOD 

Here, we give some additional details about our numerical procedures. Before performing the 
discretization of the dimensionless Eqs. (1271 ), (1301 ), (1331) and (1361 ) we choose to place ourselves 
in the frame moving at the velocity of the fluid at the outlet. Using dimensionless units, the 
velocity field in this frame is v so that in the laboratory frame u = y + v, (one should note that for 
a planar interface v = 0). Moreover, we have chosen to solve the pressure equation considering 
a perturbation of the average pressure gradient. This amounts to preconditionning the Poisson 
operator which is ill-conditionned because of the presence of high viscosity contrasts. 

The spatial discretization of equations (1271 ), (1301 ), (1351) and (1361 ) is based on finite differences 
on a staggered grid. The pressure and the phase field are evaluated at the mesh nodes, whereas 
the horizontal velocity component are evaluated at the mid-points of the horizontal links, and the 
vertical component at the mid-points of the vertical links. This staggering allows for a scheme that 
exactly guarantees mass conservation (V-v = 0). Time-stepping appears only in the phase field 
evolution equation. An explicit formulation of the time derivatives is retained, which results in the 
following sequence (the time step is indicated as a superscript): 

1. Given <f) n and p n_1 , we calculate /t™ ff (0™ , £> n_1 ) and K n . 

2. The pressure field p n is given by the solution of the equation 

For a stationnary state, the effective viscosity and the pressure are the solutions of a fixed 
point problem which converges in practice (for small enough time steps and suitable physical 
parameters). 

3. The velocity is then directly evaluated by 

r = -4- {vp n -Bn n + - m . (B2) 

4. The phase field is timestepped, 

0™+! = (f) n + dt {-W 1 ■ V0 n + e~ 2 {<f) - 3 ) n + V 2 0" + « n |V0"|} (B3) 
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Here, /2 e g is the dimensionless interpolated viscosity averaged over the gap in the most general 
case of appendix lAl Typically in our simulations, for both Newtonian and non-Newtonian fluids 
the time step was of the order 1CT 5 . 

The pressure is obtained by solving the linear system resulting from the spatial discretization 
of the Poisson equation and associated boundary conditions using The Gauss-Seidel SOR method. 
The SOR solver is initialized with the pressure field of the preceding time step, which helps to 
reduce significantly the amount of iterations needed to achieve convergence of the pressure field 
after the few initial time steps. In all our simulations, the relaxation parameter uj is set to 1,83. 
This value is chosen by trial and error in order to limit the amount of overall iterations at each 
time step and to allow for fast enough parametric studies. The convergence criterion e sor was 
chosen so that the residual of the linear system was very close Il42ll to 10~ 5 . A drawback of 
the SOR method is its sensitivity to the number of mesh points compared to more sophisticated 
methods such as, for instance, ADI. Computations performed on a 2D test-case showed that the 
number of iterations required to solve the problem on a square mesh with an imposed accuracy 
of 1CT 5 increases roughly as the number of nodes in one direction when an relaxed ADI method 
necessitates a fixed number of iterations once the problem is well resolved spatially. However, 
for the grid resolution used in our parametric studies, N x x N y — 100 x 500, the number of 
iterations varies between 20 and 50 after the few initial time steps. Not surprisingly, the number 
of iterations is found to depend on the choice of the physical parameters: more strongly shear- 
thinning fluids (higher n in Carreau's law) and higher We require more iterations to converge. I is 
also found that varying the viscosity ratio can change the number of iterations needed to converge. 
Unexpectedly, when v tends to zero the number of iterations per time step decreases, whereas 
the problem becomes numerically more challenging. This does not indicate that the SOR method 
works better for harder problems but that instead of solving the pressure field in front of the finger 
the iterations are used to solve the field in the finger which bears almost no information. The 
number of digits available in the numerical solution for the diffusion field is greatly influenced by 
the value of v. Since the overall available information in the numerical solution is set by e sor we 
adapted its value to roughly maintain constant the number of SOR iterations. In the Newtonian 
case, we found that to keep the SOR iterations value approximatively to 20 for v values of 0.05 
and 0.005, e sor needs to be set to 10~ 4 and 10 -5 , respectively. 

We conclude by a few remarks on the discretization. In order to statisfy the incompressibility 
constraint, the Poisson equation should be carefully discretized. Since the divergence operator is 
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applied to equation (1271) . all the terms of this equation should be evaluated in the middle of the 
links surrounding the point where the discrete divergence applies (see fig©. This amounts to: 



+ 1 - h (B4) 



with 



(B5) 



for the left link, the discretization associated with the other links being staightforward. Similarly 



/2ec / V ^ 2 ' J Aeff 2 'V Ax 

+ v (B6) 



with, 



and 



In order to minimize rounding errors and ensure mass conservation, the spatial discretization 
of the velocity equation must be consistent with that of the Poisson equation for the pressure. 
Staggering results in shifting the velocity components in the direction to which they relate. This 
gives: 



^=- ( S^- +b O^ <B9) 

V U - + <r>^ + 5 U=b (B10) 



The treatment of the phase field equation is straightforward and analogous to that found in 
Ref. [21] except for the advective term. The velocity on the the nodes is recovered by a linear 
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interpolation: 



2 V *+3>J l -2.-? y Ax 
^(^■^.J ^ 1 ^- 1 . (BID 

The use of the phase field to compute a continuous equivalent of the curvature k requires to intro- 
duce a numerical cutoff to avoid infinite values at centers of curvature. Wherever | V0| is smaller 
than 1CT 4 , k is replaced by 0. 
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We = 0.005 




Figure 4: Maps of the viscosity for the two-plateau law with a = 0.3 at various Weissenberg numbers. 
Darker tones correspond to more viscous regions. 
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Figure 5: Finger width A versus We at fixed values of T and a for the two-plateau viscosity law. 
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Figure 6: Selected finger width as a function of B for the effective viscosity of Eq. (I4TT) with three different 
values of a, and comparison to the Newtonian case; for all simulations, v = 5 X 10~ 3 , e=0.02, Ax = 0.01. 
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Figure 7: Finger width versus We for the one-plateau law, either at fixed T, or at fixed T out , which is the 
dimensionless surface tension computed using the effective viscosity at the outlet [see Eq. (f44ll. 




~i 1 r 



i r 



J* 



A' 



_L 



w 4cm b .66mm 

w 4cm b .25mm 

w 8cm b .5mm 

We 20 

We 100 
i i i_ 



▲ 

-B- 



..kfl 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 

r/2(We) 1/2 



Figure 8: Finger width as a function of the reduced parameter We n r for two sets of simulations, and 
comparison to the experimental data. 
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Figure 9: Staggering of the fields 
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